######################################################################
#Replication file for "The Impact of China's AIIB on the World Bank"
#Authors: Jing Qian, James Vreeland, Jianzhi Zhao
#Codes to replicate Figure A.5 - Figure A.7
######################################################################
rm(list = ls())

#To install the bpCausal package (ver 0.0.1)
# install.packages('devtools', repos = 'http://cran.us.r-project.org') # if not already installed
# devtools::install_github('liulch/bpCausal')

#Load bpCausal package
library(bpCausal)

#Other packages used
library(tidyverse)
library(RColorBrewer)

#Setwd
setwd("C:/Users/qianj/Dropbox (Princeton)/AIIB_WB_Replication")

#Load custom functions
#Note: Given the number of tests performed in the appendix, 
#several custom functions are created to implement estimation methods in a more efficient way.
#These functions merely create a wrap-up of existing functions, like bpCausal and gsynth,
#in order to avoid typing parameters repeatedly.

source("code/custom_functions.R")

###################################
#Results for Figure A.5-A.7
###################################
#----------------
#(1) Load original results
#----------------
load("result/bpcausal_original.RData")


###################################
#Produce Figure A.5
###################################
#------------------
#(1) Prepare results
#------------------
#Groups
iso2c.tr = result.original$raw.id.tr %>%
  unique() #All treated id
iso2c.nochina = setdiff(iso2c.tr, "CN") #Exclude China
iso2c.nobrza = setdiff(iso2c.tr, c("BR", "ZA")) #Exclude Brazil and South Africa

#Summarize
result.use = effSummary(result.original)
result.nochina.use = effSummary(result.original,
                                usr.id = iso2c.nochina)
result.nobrza.use = effSummary(result.original,
                               usr.id = iso2c.nobrza)

#------------------
#(2) Parameters
#------------------
n = 3
xlim = c(-2, 2)
cex = 2
xtext = 0.6
colors = brewer.pal(n = n, 
                    name = "Dark2")

#--------------------
#(3) Plot
#--------------------
#----------
#Setup
png(filename = "figure/Figure_A5.png",
    height = 800,
    width = 1200)

par(mar = c(5, 1, 5, 1))

#----------
#Empty plot
plot(1,
     type = "n",
     xlab = "Estimated Average Treatment Effect on the Treated",
     ylab = "",
     yaxt = "n",
     xaxt = "n",
     xlim = xlim,
     ylim = c(n + 3 , 0),
     cex.lab = 2)
box()

#-------------
#Original
#-------------
points(x = result.use$est.avg[1],
       y = 2,
       col = colors[1],
       pch = 19,
       cex = cex)
arrows(x0 = result.use$est.avg[2],
       x1 = result.use$est.avg[3],
       y0 = 2,
       y1 = 2,
       length = 0.05,
       angle = 90,
       code = 3,
       col = colors[1],
       cex = cex,
       lwd = cex)

#----------------
#Remove outliers
#----------------
#Remove China
points(x = result.nochina.use$est.avg[1],
       y = 2 + 1,
       col = colors[1+1],
       pch = 19,
       cex = cex)
arrows(x0 = result.nochina.use$est.avg[2],
       x1 = result.nochina.use$est.avg[3],
       y0 = 2 + 1,
       y1 = 2 + 1,
       length = 0.05,
       angle = 90,
       code = 3,
       col = colors[1+1],
       cex = cex,
       lwd = cex)

#Remove Brazil & South Africa
points(x = result.nobrza.use$est.avg[1],
       y = 2 + 2,
       col = colors[1+2],
       pch = 19,
       cex = cex)
arrows(x0 = result.nobrza.use$est.avg[2],
       x1 = result.nobrza.use$est.avg[3],
       y0 = 2 + 2,
       y1 = 2 + 2,
       length = 0.05,
       angle = 90,
       code = 3,
       col = colors[1+2],
       cex = cex,
       lwd = cex)

#---------------
#Other elements
#Abline
abline(v = 0,
       lty = "dashed",
       lwd = cex,
       col = "grey")
#Axis
axis(side = 1,
     at = seq(min(xlim),
              max(xlim),
              1),
     cex.axis = cex)

#Title
title("World Bank Infrastructure Projects",
      cex.main = cex,
      line = 1)

#Legend
legend("topright",
       legend = c("Original",
                  "Exclude China",
                  "Exclude BR & ZA"),
       cex = cex,
       col = colors[1:3],
       lty = rep(1, n),
       bty = "n",
       lwd = cex)

#-----------------
#Close
dev.off()


###################################
#Produce Figure A.6
###################################
#------------------
#(1) Prepare results
#------------------
#Groups
iso2c.tr = result.original$raw.id.tr %>%
  unique() #All treated id

#Summarize
result.jk.use = list() #container
for (i in 1:length(iso2c.tr)){
  result.jk.use[[i]] = effSummary(result.original,
                                  usr.id = iso2c.tr[-i])
}
names(result.jk.use) = iso2c.tr

#Data
est.use = result.jk.use %>%
  lapply(., function(x) x$est.avg) #Average ATT

coef.use = est.use %>%
  lapply(., function(x) x[1]) %>%
  unlist() #Coefficient

ci.l.use = est.use %>%
  lapply(., function(x) x[2]) %>%
  unlist() #lower bound

ci.u.use = est.use %>%
  lapply(., function(x) x[3]) %>%
  unlist() #upper bound

#-------------
#(2) Parameters
#-------------
n = length(iso2c.tr)
xlim = c(-2, 2)
cex = 2
xtext = 0.6
colors = brewer.pal(n = 3,
                    name = "Dark2")

#---------------
#(3) Plot
#---------------
#Setup
png(filename = "figure/Figure_A6.png",
    width = 800,
    height = 1200)

par(mar = c(5, 5, 5, 1))

#Empty plot
plot(1, 
     type = "n",
     xlab = "Estimated Average Treatment Effect on the Treated",
     ylab = "",
     yaxt = "n",
     xaxt = "n",
     xlim = xlim,
     ylim = c(n + 0.5, 0.5),
     cex.lab = 2)
title(main = "World Bank Infrastructure Projects",
      line = 1,
      cex.main = cex)
box()

#Points
points(x = coef.use,
       y = 1:n,
       col = colors[2],
       pch = 19,
       cex = cex)
#Arrows
arrows(x0 = ci.l.use,
       x1 = ci.u.use,
       y0 = 1:n,
       y1 = 1:n,
       length = 0.05,
       angle = 90,
       code = 3,
       col = colors[2],
       cex = cex,
       lwd = cex)

#Other elements
#abline
abline(v = 0,
       lty = "dashed",
       lwd = cex,
       col = "gray")
#Axis
axis(side = 1,
     at = seq(min(xlim),
              max(xlim),
              1),
     cex.axis = cex)
#Removed countries
axis(side = 2,
     at = 1:n,
     labels = names(result.jk.use),
     xpd = NA,
     cex.axis = 2,
     las = 2,
     tck = 0.000001)
#-----------
#Close
dev.off()

